##Get data ready for the recursive example
#

rm(list=ls())
library(ctgt)
library(stringr)
library(xlsx)


n= 100
m = 5
X = matrix(0, n, m,byrow = T )
for ( i in 1:n){
  set.seed(1234+i)
  X[i,] =  as.vector(arima.sim(model = list(order = c(1, 0, 0), ar = 0.2), n = m) )
}

y = rbinom(n,1,0.6)
X[which(y==1),1:3] = X[which(y==1),1:3] + 0.8


#cor(X[,1:50],X[,1:50])
library(stringr)
xs = str_replace_all(paste(rep("x",m),seq(1,m,1)),fixed(" "), "")
colnames(X) = xs

#Rprof()
actgt(y=y,X=X,xs=xs,hyps=c("x3"))
#Rprof(NULL)
#summaryRprof()

sqrW = sqrt(mean(y)*(1-mean(y)) ) # sqrt of covariance of y

WIHZ = sqrW *(sweep(X,2,colMeans(X))) ## W^{1/2}*(I-H)Z
IHZ = WIHZ/sqrW ##(I-H)Z
# the full model
Tf = sum(colSums(y*IHZ[,xs,drop=F])^2)# full test statistic = y^t (I-H) Z Z^t (I-H) y
Lamf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # eigen(W^{1/2}*(I-H)Z%*%Z^t%*%(I-H)%*%W^{1/2}) 
Cf = criticalvalue(Lamf )

Tf

Cf

tt = sapply(1:m, function(x) sum(colSums(y*IHZ[,x,drop=F])^2) )

ww = sapply(1:m, function(x) sum(round(eigen(tcrossprod(WIHZ[,x,drop=F]),symmetric = T,only.values = T)$values,8) ) )

tt/ww

namebank = c("3", "13", "23", "34", "35","123", "134", "135", "234", "235", "345", "1234", "1235", "1345", "2345", "12345")
sublis = list(3, c(1,3), c(2,3), c(3,4),c(3,5), c(1,2,3), c(1,3,4),c(1,3,5), c(2,3,4),c(2,3,5),c(3,4,5), c(1,2,3,4), c(1,2,3,5), c(1,3,4,5),c(2,3,4,5),c(1,2,3,4,5))
tmin = c()
level = c()
for(i in 1:16){
  tmin[i] = sum(colSums(y*IHZ[,sublis[[i]],drop=F])^2)#
  level[i] = sum(round(eigen(tcrossprod(WIHZ[,sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) )
}

## cmax 
laf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # 
las = round(eigen(tcrossprod(WIHZ[,3,drop=F]),symmetric = T,only.values = T)$values,8) # 
levs = seq(sum(las), sum(laf), length.out = 200)
cmax=c()
for (i in 1:length(levs)){
  cmax[i] = criticalvalue(getL(laf,las,levs[i]))
}
cmax

## cture
#sublis = list(3, c(1,3), c(2,3), c(3,4),c(3,5), c(1,2,3), c(1,3,4),c(1,3,5), c(2,3,4),
#c(2,3,5),c(3,4,5), c(1,2,3,4), c(1,2,3,5), c(1,3,4,5),c(2,3,4,5),c(1,2,3,4,5))
crt = c()
for (i in 1: 16){
  la = round(eigen(tcrossprod(WIHZ[, sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) 
  crt[i] = criticalvalue(la)
}
tc3 = list(tmin = tmin,level = level,cmax =cmax, crt = crt,levs = levs)
save(tc3, file="gmincmax3.RData")

########


########

rm(list=ls())

n= 100
m = 5
X = matrix(0, n, m,byrow = T )
for ( i in 1:n){
  set.seed(1234+i)
  X[i,] =  as.vector(arima.sim(model = list(order = c(1, 0, 0), ar = 0.2), n = m) )
}

y = rbinom(n,1,0.6)
X[which(y==1),1:3] = X[which(y==1),1:3] + 0.8


#cor(X[,1:50],X[,1:50])
library(stringr)
xs = str_replace_all(paste(rep("x",m),seq(1,m,1)),fixed(" "), "")
colnames(X) = xs

actgt(y=y,X=X,xs=xs,hyps=c("x2"))



sqrW = sqrt(mean(y)*(1-mean(y)) ) # sqrt of covariance of y

WIHZ = sqrW *(sweep(X,2,colMeans(X))) ## W^{1/2}*(I-H)Z
IHZ = WIHZ/sqrW ##(I-H)Z
# the full model
Tf = sum(colSums(y*IHZ[,xs,drop=F])^2)# full test statistic = y^t (I-H) Z Z^t (I-H) y
Lamf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # eigen(W^{1/2}*(I-H)Z%*%Z^t%*%(I-H)%*%W^{1/2}) 
Cf = criticalvalue(Lamf )

Tf

Cf

tt = sapply(1:m, function(x) sum(colSums(y*IHZ[,x,drop=F])^2) )

ww = sapply(1:m, function(x) sum(round(eigen(tcrossprod(WIHZ[,x,drop=F]),symmetric = T,only.values = T)$values,8) ) )

tt/ww

namebank = c("2", "12", "23", "24", "25","123", "124", "125", "234", "235", "245", "1234", "1235", "1245", "2345", "12345")
sublis = list(2, c(1,2), c(2,3), c(2,4),c(2,5), c(1,2,3), c(1,2,4),c(1,2,5), c(2,3,4),c(2,3,5),c(2,4,5), c(1,2,3,4), c(1,2,3,5), c(1,2,4,5),c(2,3,4,5),c(1,2,3,4,5))
tmin = c()
level = c()
for(i in 1:16){
  tmin[i] = sum(colSums(y*IHZ[,sublis[[i]],drop=F])^2)#
  level[i] = sum(round(eigen(tcrossprod(WIHZ[,sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) )
}

## cmax 
laf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # 
las = round(eigen(tcrossprod(WIHZ[,3,drop=F]),symmetric = T,only.values = T)$values,8) # 
levs = seq(sum(las), sum(laf), length.out = 200)
cmax=c()
for (i in 1:length(levs)){
  cmax[i] = criticalvalue(getL(laf,las,levs[i]))
}
cmax

## cture
#sublis = list(3, c(1,3), c(2,3), c(3,4),c(3,5), c(1,2,3), c(1,3,4),c(1,3,5), c(2,3,4),
#c(2,3,5),c(3,4,5), c(1,2,3,4), c(1,2,3,5), c(1,3,4,5),c(2,3,4,5),c(1,2,3,4,5))
crt = c()
for (i in 1: 16){
  la = round(eigen(tcrossprod(WIHZ[, sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) 
  crt[i] = criticalvalue(la)
}

tc2 = list(tmin = tmin,level = level,cmax =cmax, crt = crt,levs = levs)
save(tc2, file="gmincmax2.RData")




########


########

rm(list=ls())


n= 100
m = 5
X = matrix(0, n, m,byrow = T )
for ( i in 1:n){
  set.seed(1234+i)
  X[i,] =  as.vector(arima.sim(model = list(order = c(1, 0, 0), ar = 0.2), n = m) )
}

y = rbinom(n,1,0.6)
X[which(y==1),1:3] = X[which(y==1),1:3] + 0.8


#cor(X[,1:50],X[,1:50])
library(stringr)
xs = str_replace_all(paste(rep("x",m),seq(1,m,1)),fixed(" "), "")
colnames(X) = xs

actgt(y=y,X=X,xs=xs,hyps=c("x1"))
actgt(y=y,X=X,xs=xs,hyps=c("x1"),maxit = 10)


sqrW = sqrt(mean(y)*(1-mean(y)) ) # sqrt of covariance of y

WIHZ = sqrW *(sweep(X,2,colMeans(X))) ## W^{1/2}*(I-H)Z
IHZ = WIHZ/sqrW ##(I-H)Z
# the full model
Tf = sum(colSums(y*IHZ[,xs,drop=F])^2)# full test statistic = y^t (I-H) Z Z^t (I-H) y
Lamf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # eigen(W^{1/2}*(I-H)Z%*%Z^t%*%(I-H)%*%W^{1/2}) 
Cf = criticalvalue(Lamf )

Tf

Cf

tt = sapply(1:m, function(x) sum(colSums(y*IHZ[,x,drop=F])^2) )

ww = sapply(1:m, function(x) sum(round(eigen(tcrossprod(WIHZ[,x,drop=F]),symmetric = T,only.values = T)$values,8) ) )

tt/ww

namebank = c("1", "12", "13", "14", "15","123", "124", "125", "134", "135", "145", "1234", "1235", "1245", "1345", "12345")
sublis = list(1, c(1,2), c(1,3), c(1,4),c(1,5), c(1,2,3), c(1,2,4),c(1,2,5), c(1,3,4),c(1,3,5),c(1,4,5), c(1,2,3,4), c(1,2,3,5), c(1,2,4,5),c(1,3,4,5),c(1,2,3,4,5))
tmin = c()
level = c()
for(i in 1:16){
  tmin[i] = sum(colSums(y*IHZ[,sublis[[i]],drop=F])^2)#
  level[i] = sum(round(eigen(tcrossprod(WIHZ[,sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) )
}


## cmax 
laf = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # 
las = round(eigen(tcrossprod(WIHZ[,3,drop=F]),symmetric = T,only.values = T)$values,8) # 
levs = seq(sum(las), sum(laf), length.out = 200)
cmax=c()
for (i in 1:length(levs)){
  cmax[i] = criticalvalue(getL(laf,las,levs[i]))
}
cmax

## cture
#sublis = list(3, c(1,3), c(2,3), c(3,4),c(3,5), c(1,2,3), c(1,3,4),c(1,3,5), c(2,3,4),
#c(2,3,5),c(3,4,5), c(1,2,3,4), c(1,2,3,5), c(1,3,4,5),c(2,3,4,5),c(1,2,3,4,5))
crt = c()
for (i in 1: 16){
  la = round(eigen(tcrossprod(WIHZ[, sublis[[i]],drop=F]),symmetric = T,only.values = T)$values,8) 
  crt[i] = criticalvalue(la)
}

tc1 = list(tmin = tmin,level = level,cmax =cmax, crt = crt,levs = levs)
save(tc1, file="gmincmax1.RData")